function [model] = ContW(data,opts)
% [model] = ContW(data,opts): Efficient Building W
% Input:
%       - data: the data matrix of size nSmp x nFea, where each row is a sample
%               point
%       opts: options for this algorithm
%           - p: the number of landmarks picked (default 1000)
%           - r: the number of nearest landmarks for representation (default 5)
%           - a: weight in manifold ranking, score = (I - aS)^(-1)y, default  0.99
%           - mode: landmark selection method, currently support
%               - 'kmeans': use centers of clusters generated by kmeans (default)
%               - 'kmeans_one_of_ten': use 10% data samples
%               - 'random': use randomly sampled points from the original
%                           data set 
%               - 'given': use given anchors in opts.ancs 
%           The following parameters are effective ONLY in mode 'kmeans'
%           - kmNumRep: the number of replicates for initial kmeans (default 1)
%           - kmMaxIter: the maximum number of iterations for initial kmeans (default 5)
%
% Output:
%       - model: the learned model 
%
%
%   Written by Bin Xu (binxu986 AT gmail.com)



% Set and parse parameters
if (~exist('opts','var'))
   opts = [];
end

p = 1000;
if isfield(opts,'p')
    p = opts.p;
end

r = 5;
if isfield(opts,'r')
   r = opts.r;
end

a = 0.99;
if isfield(opts,'a')
   a = opts.a;
end

mode = 'kmeans';
if isfield(opts,'mode')
    mode = opts.mode;
end
kmMaxIter = 5;
if isfield(opts,'kmMaxIter')
    kmMaxIter = opts.kmMaxIter;
end
kmNumRep = 1;
if isfield(opts,'kmNumRep')
    kmNumRep = opts.kmNumRep;
end
nSmp =size(data,1);

% Landmark selection
if strcmp(mode,'kmeans')
    [~,landmarks]=litekmeans(data,p,'MaxIter',kmMaxIter,'Replicates',kmNumRep);
elseif strcmp(mode, 'kmeans_one_of_ten')
    indSmp = randperm(nSmp);
    [~,landmarks]=litekmeans(data(indSmp(1:ceil(nSmp/10)),:),p,'MaxIter',kmMaxIter,'Replicates',kmNumRep);
    clear indSmp;
elseif strcmp(mode,'random')
    indSmp = randperm(nSmp);
    landmarks = data(indSmp(1:p),:);
    clear indSmp
elseif strcmp(mode,'given')
    landmarks = opts.ancs;
else
    error('mode does not support!');
end

model.p=p;
model.landmarks = landmarks;
model.a = a;
model.r = r;
model.mode = mode;

%%% Z construction
%t0 = clock;
Z = [];
NumOneIter = 1e5;
for j = 1:floor(nSmp/NumOneIter)
    fea= data(1:NumOneIter, :  );
    D = EuDist2(fea,landmarks);
    ItSmp = NumOneIter;
    dump = zeros(ItSmp,r);
    idx = dump;
    for i = 1:r
        [dump(:,i),idx(:,i)] = min(D,[],2);
        temp = (idx(:,i)-1)*ItSmp+[1:ItSmp]';
        D(temp) = 1e10;
    end
    dump = bsxfun(@rdivide,dump,dump(:,r));
    dump = 0.75 * (1 - dump.^2);
    Gsdx = dump;
    Gidx = repmat([1:ItSmp]',1,r);
    Gjdx = idx;
    Z=[ Z ;sparse(Gidx(:),Gjdx(:),Gsdx(:),ItSmp,p)];
    data(1:NumOneIter, :  ) = [];
end
if(~isempty(data))
    D = EuDist2(data, landmarks);
    ItSmp = size(data,1);
    dump = zeros(ItSmp,r);
    idx = dump;
    for i = 1:r
        [dump(:,i),idx(:,i)] = min(D,[],2);
        temp = (idx(:,i)-1)*ItSmp+[1:ItSmp]';
        D(temp) = 1e10;
    end
    dump = bsxfun(@rdivide,dump,dump(:,r));
    dump = 0.75 * (1 - dump.^2);
    Gsdx = dump;
    Gidx = repmat([1:ItSmp]',1,r);
    Gjdx = idx;
    Z=[ Z ;sparse(Gidx(:),Gjdx(:),Gsdx(:),ItSmp,p)];
    % t1 = clock; time1 = etime(t1, t0)           %time  <---------------
    clear dump;
end
% Learning the model for W, W = H*H';
feaSum = (sum(Z,1));
D = Z*feaSum';
D = max(D, 1e-12);
D = D.^(-.5);
H = spdiags(D,0,nSmp,nSmp)*Z;
model.feaSum = feaSum;
model.H = H;



